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ABSTRACT 

Detection of -B-mode polarization of the cosmic microwave background (CMB) radiation is one of the 
frontiers of observational cosmology. Because they are an order of magnitude fainter than _E-modes, it 
is quite a challenge to detect 73-modes. Having more manageable systematics, interferometers prove 
to have a substantial advantage over imagers in detecting such faint signals. Here, we present a 
method for Bayesian inference of power spectra and signal reconstruction from interferometric data of 
the CMB polarization signal by using the technique of Gibbs sampling. We demonstrate the validity 
of the method in the flat-sky approximation for a simulation of an interferometric observation on a 
finite patch with incomplete uw-plane coverage, a finite beam size and a realistic noise model. With 
a computational complexity of 0(n 3 / 2 ), n being the data size, Gibbs sampling provides an efficient 
method for analyzing upcoming cosmology observations. 

Subject headings: cosmology:observations, cosmic microwave background, polarization 
instrumentation:intcrfcromctcrs, methods: data analysis, methods: statistical 



1. INTRODUCTION 

The cosmic microwave background (CMB) polariza- 
tion signal can be decomposed into a scalar E compo- 
nent and a pseudo-scalar B component (Zaldarriaga & 
Seljak 1997; Kamionkowski et al. 1997). The largest 
contribution to the CMB polarization comes from the 
scalar metric perturbations produced by density fluctua- 
tions, which produce only i_-type polarization. At small 
angular scales (£ <~ 1000) gravitational lensing due to 
large-scale structure transforms a small portion of the E- 
modes into B-modes (Zaldarriaga & Seljak 1998). The 
more interesting source of _3-type polarization is the pri- 
mordial tensor metric perturbations produced by grav- 
itational waves created during inflation (Zaldarriaga & 
Seljak 1997; Kamionkowski et al. 1997). Since tensor 
modes dominate on large angular scales (I ~ 100), detec- 
tion of B-modes at these scales offers an excellent probe 
for the inflationary epoch whose energy scale is propor- 
tional to the amplitude of primordial gravitational waves 
(Hu & White 1997). 

Because J5-modes are not produced by scalar pertur- 
bations, they are smaller than E- modes by more than 
an order of magnitude. Detection of such weak signals, 
at a level of tensor-to-scalar ratio of 0.01, requires ex- 
cellent control of systematic effects. Since traditional 
imagers measure Q and U Stokes parameters by differ- 
encing two orthogonal polarizations, mismatched beams 
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and pointing errors cause leakage from the much stronger 
temperature signal into the Q and U signals, signifi- 
cantly contaminating the much weaker B-modes (Hu et 
al. 2003). Interferometers, on the other hand, directly 
measure the Stokes parameters without subtraction of 
the signals from different detectors. Thus, mismatch 
in the beam patterns or differential pointing errors do 
not cause contamination of the polarization by the tem- 
perature signal (Bunn 2007). Morover, for finite sky 
patches and pixellated maps, E and i3-modes mix into 
each other, causing major contamination of -B-modes by 
much stronger E- modes (Lewis et al. 2002; Bunn 2003). 
Since interferometric data live in Fourier space, separa- 
tion of E and 73-modes can be achieved more cleanly by 
interferometers than imagers (Park et al. 2003; Park & 
Ng 2004). 

Interferometers have already been applied to the de- 
tection of the polarized CMB signal. The first detection 
of polarization anisotropies in the CMB was achieved by 
DASI (Kovac et al. 2002). CBI (Pearson et al. 2003) 
and VSA (Dickinson et al. 2004; Grainge et al. 2003) 
obtained detailed observations enabling them to extract 
the S-mode polarization angular power spectrum up to 
I ~ 600. Advancing techniques, such as bolometric inter- 
ferometry employed by the QUBIC experiment (Battis- 
telli et al. 2010), provide promising developments in de- 
tecting the long-sought B-mode polarization anisotropies 
yet to be observed. 

In comparison to alternative methods of extracting 
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power spectra such as maximum likelihood and pseudo- 
Ci estimators, which often scale as 0(n 3 ) and 0(n 3 / 2 ) 
respectively, the method of Gibbs sampling (Jewell et al. 
2004; Wandelt et al. 2004) has an advantage in dealing 
with the demands of future cosmology observations be- 
cause it provides simultaneous estimation of power spec- 
trum and signal by sampling them from the joint poste- 
rior probability density through a Markov chain Monte 
Carlo process with 0(n 3 / 2 ) computational complexity. 
Gibbs sampling has already been used to analyze the 
WMAP temperature (O'Dwyer et al. 2004; Dickinson 
et al. 2009; Larson et al. 2011) and polarization (Lar- 
son et al. 2007; Eriksen et al. 2007; Komatsu et al. 
2011) data. Sutter et al. (2011) examined the applica- 
tion of Gibbs sampling to interferometric observations of 
the CMB temperature signal. 

Here, we will investigate the application of Gibbs sam- 
pling to the polarized CMB signal observed by an inter- 
ferometer. Our analysis is an extension of Gibbs sam- 
pling as applied to interferometric data (Sutter et al. 
2011) to polarized signals along the lines of Larson et 
al. (2007). In Section 2 we investigate the method of 
Gibbs sampling as applied to interferometric polarime- 
try. In Section 3. we discuss the simulation of interfer- 
ometric observation of the polarized CMB signal on a 
finite patch in the flat-sky approximation. In Section 4 
we present our results of polarized power spectra and sig- 
nal reconstructions. Finally, in Section 5 we make some 
comments and concluding remarks. 

2. METHOD OF GIBBS SAMPLING 

In the flat-sky approximation, we describe the CMB 
signal dimensional vector, s, of the Fourier 

transform of the discretized sky maps of n p pixels; 

s = (. . . ,Ti,Ei,Bi, . ..); i = 0, ...,n p — 1, where / de- 
notes the Fourier transform of /. The covariance matrix 
S = (s s t) of the CMB signal is a block diagonal matrix 
with a 3 x 3 submatrix Cj at each pixel i: 
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where li = 27r|{ti| and itj is the position vector of the i th 
pixel in the Fourier plane. Larson et al. (2007) have ap- 
plied Gibbs sampling to the full-sky WMAP polarization 
data. Their analysis can be extended to interferometric 

observations by describing the visibility data, d, from a 
polarimetric observation as pixelated maps of the Stokes 
parameters T, Q, and U : 

d = H R s + n (2) 

where n is a Gaussian realization of the noise, H is a 
linear operator that includes convolution with an instru- 
ment beam A and R, a block diagonal matrix with a 
3x3 submatrix R^ at each pixel, is the transformation 
of the T, E, and B components of the signal s into the 
Fourier transform of the Stokes parameters. For an in- 
terferometric data set, in the flat-sky approximation, H 
and R, can be written as 
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where F is a Fourier transform operator, / is an inter- 
ferometer patten in the uw-plane, and <pi is the angular 
position of the i th pixel in the uu-plane. 

Larson et al. (2007) investigated Gibbs sampling as ap- 
plied to polarized signals. The principle approach is to 
sample s and S from the joint density P(S, s, d ), which 
can be obtained by a Markov chain Monte Carlo process 
by successively sampling from the conditional distribu- 
tions P(s | S, d) and P(S | a, d) oc P(S | s ). Starting 
from an initial guess S°, sampling is done in an iterative 
fashion by (Larson et al. 2007) 
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After some "burn-in" steps the stationary distribution 
of the Markov chain is reached and the samples approx- 
imate to being samples from the sought-after joint dis- 
tribution. Sampling from the joint distribution by this 
technique is called the method of Gibbs sampling. 
Given the current covariance matrix S a , the sky signal, 
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y, separated into the mean field, x, and 



fluctuation, y, parts, is sampled by solving the following 
equations (Larson et al. 2007): 
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(R S" 1 R T + H T N" 1 H)y = R S~ 1/2 f+H T r\T 1/2 x (8) 

where N = (n n ^) is the noise covariance matrix, which 
is diagonal (White et al. 1999) and has entries equal to 
Nij = vfSij where ^ is the noise variance for the i th 

pixel in the ww-plane, and £ and \ are Gaussian random 
maps having values of zero mean and unit variance in 
each pixel for each of the T, Q, and U components. 

We obtain numerical solutions for the Eq. [7] and Eq. [8] 
by the preconditioned conjugate-gradient method (Press 
et al. 1986). A good choice for the preconditioner is the 
inverse of the diagonal part of the operator: 
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The noise portion of the preconditioner, Pat, can be writ- 
ten as (Sutter et al. 2011) 
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where A is the Fourier transform of the beam pattern. 

The signal polarization map, s a , sampled from 
P( s | S a_1 , d ), is used to sample the signal covariance 
matrix from P{ S | s a ) by computing the unnormalized 
variance at in an annulus of radius £/2ir. We can de- 
fine uniform bins b = [£ m in, £max] in which Ctl{t + 1) is 
roughly constant. Then ag is defined for bin b as (Larson 
et al. 2007) 
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where Si = (Ti,Ei,Bi) is a three- vector at the i th pixel. 
Sampling from the probability density P ( | s a ), 
which is an inverse Wishart distribution with mb degrees 
of freedom, can be done by drawing m& = Pb — 2 (as- 
suming a Jeffreys' ignorance prior) vectors from a Gaus- 
sian distribution with covariance matrix , where pb 
is the number of pixels in the bin b. The required sam- 
ple Cb is, then, the inverse of the sum of outer products 
of these independently sampled vectors (Larson et al. 
2007). The actual power spectrum coefficients are given 

by c t = c b /£(e+i). 

Following Sutter et al. (2011), the Gelman-Rubin (G- 
R) statistic is employed to determine that the stationary 
distribution of the Markov chain has been reached. Given 
multiple instances of chains, convergence is reached when 
the potential scale reduction factor of the G-R statistic, 
determined by the ratio of the variance within each chain 
to the variance among chains, assumes a value less than 
a given tolerance for each bin (Gelman & Rubin 1992). 

3. SIMULATIONS 

We construct the input Q and U maps by transform- 
ing realizations of E and B signals over 10-degree square 
patches with 128 pixels per side. The realizations are 
created as maps of Gaussian fluctuations with a covari- 
ance Cj (Eq. 1) at each pixel whose components are 
produced by CAMB (Lewis et al. 2000). The cosmo- 
logical parameters used for CAMB are consistent with 
the 7-year results of WMAP (Larson et al. 2011; Ko- 
matsu et al. 2011); fi M = 0.27, fl A = 0.73, Cl b = 0.045, 
and H — TQkms" 1 Mpc -1 . The tensor-to-scalar ratio is 
taken to be T/S = 0.01. The angular resolution of the 
signal maps is 4.7 arcminutes corresponding to a max- 
imum available multipole of l m ax — 2 300. The spatial 
resolution in the ui>-plane is 5.73 A. Although the patch 
size is too large to employ the fiat-sky approximation, it 
is still useful in exploring the validity of our technique. 

The primary beam pattern A is modeled as a Gaus- 
sian with peak value of unity and standard deviation of 
2.5 degrees allowing us to include all Fourier modes up 
to the Nyquist frequency in the analysis. Although a 
smaller beam size would further reduce the edge-effects 
caused by the periodic boundary conditions of the fast 
Fourier transformations, it would also require a longer 
computation time. 

The interferometer is constructed by randomly placing 
16 antennas with diameters of 10 cm in the ww-plane and 
uniformly rotating the baselines over a period of 12 hours 
while observing the same sky patch. The observation 
frequency is 30 GHz with a 10-GHz bandwidth. With 
this frequency and antenna radius the minimum avail- 
able multipole is £ m i n — 28. The interferometer pattern 
/ is constructed by placing a value of one at each pixel 
that coincides with a baseline length during the obser- 
vation period and zeros everywhere else. With the given 
number of antennas and pixel resolution, the resulting 
interferometer pattern provides us with a fairly realistic 
case of incomplete uw-plane coverage as shown in Fig- 
ure 1. With this configuration, the ww-plane coverage is 
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Figure 1. Interferometer pattern, I, created over an observation 
period of 12 hours by 16 antennas of radius 5A randomly placed in 
the mi-plane. 
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Figure 2. Coverage of itu-plane for each bin A;, = [l™ m , £™ ax ]. 
Shown are the percentages of pixels intersected by baseline vec- 
tors during the 12-hour observation period in each bin. A pixel 
in Fourier space is said to be in bin b if its position vector u p i x 
satisfies 2n\u p i x \ £ A5. 



roughly 70%, which varies for each £-bin, as shown in 
Figure 2. 

The noise at each pixel for the temperature data is 
obtained from the total observation time that all base- 
lines spend in the pixel. The noise variance is given as 
1 / \J t° hs . The overall temperature signal-to-noise 
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ratio is set to 50 by scaling the noise variance at each 
pixel by the constant |H R s|/50|n|. Gaussian realiza- 
tions of this noise are used to create the T, Q, and U 
data. The corresponding signal-to-noise ratios for Q and 
U signals become 2.2 and 2.3 respectively. 

Our analysis had four independent chains, each chain 
having 2 200 iterations of which the first 200 are dis- 
carded for a "burn-in" phase. The G-R statistic reached 
less than 1.20 for each bin in about 83 hours. The com- 

1 Gelman et al. (2004) suggests that values below 1.2 are "ac- 
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Figure 3. The percentage of uv-plane coverage versus the number 
of steps, after the burn-in phase, required to reach convergence at 
each £-b'm. 



putation spent 150MB of memory on 4 cores of an Intel 
dual six-core X5650 Westmere 2.66GHz machine. 

Figure 3 shows the wu-plane coverage versus the num- 
ber of steps, after burn-in, required to reach convergence 
for each bin. We see that convergence time and uv- 
plane coverage are weakly correlated. Incomplete cov- 
erage leads to a larger correlation length for small power 
leading to a longer convergence time. Therefore, low- 
coverage bins have larger effect on overall performance, 
as expected. 

4. RESULTS 
4.1. Power Spectra 

The mean posterior power spectra of the four indepen- 
dent chains, together with the associated uncertainties 
at each €-bin obtained after convergence is reached, is 
shown in Figure 4. The input power spectra, which are 
used to construct signal realization, and the spectra of 
the signal realization are also shown in Figure 4. Nearly 
all of our estimates fall within la of the expected value, 
and none of them is outside of the 2a width. 

Although the effect of the incomplete wu-plane cover- 
age is not evident in Figure 4, if we consider the sizes of 
the uncertainties, relative to the mean posterior, we see 
that at the bins with weak ww-plane coverage the rela- 
tive sizes of the uncertainties are larger whereas at the 
bins with complete coverage the relative sizes are smaller 
as expected. The effect of sample variance dominates at 
low £ values where the sizes of the uncertainties, relative 
to the mean posterior, are larger due to the finite size of 
the sky patch. 

Gibbs sampling also provides higher-order statistical 
information such as the two-point correlations between 
^-bins. Off-diagonal components of the correlation ma- 
trices for TE 1 EE and BB power spectra are shown in 
Figure 5. There is a slight correlation between adjacent 
bins, which is the result of reduced Fourier space reso- 
lution caused by finite beam width. The correlation is 
more pronounced at high I and low signal-to-noise ra- 
tio, as seen in BB correlation matrix of Figure 5. Since 
the correlation matrices carry information about data 
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regions larger than bin sizes, the power spectra are over- 
sampled causing anti-correlation between nearby bins. 
Since the power in a region is constraint by the data, 
whenever a large value is sampled at a certain bin values 
of samples from the other bins in the same data region 
are reduced (Eisner & Wandelt 2012). 

Samples of power spectra produced by Gibbs sampling 
have highly non-Gaussian probability densities. As an 
example, we show marginalized posterior joint distribu- 
tions of EE and TE power spectra for various £-bins in 
Figure 6. Although combining many modes into each 
bin has an overall Gaussianizing effect, non-Gaussianity 
of the distributions is still clearly visible. 

4.2. Signal Reconstruction 

We compute Wiener-filtered maps by solving Eq. 7 
for x a . The mean value of these maps (F _1 Rx), trans- 
formed into Stokes variables T, Q and U and averaged 
over all iterations, are shown in Figure 7. The Wiener 
filter provides the information content of the data by fil- 
tering out the fluctuations outside of the primary beam, 
where there is no data. The fluctuations obtained by 
solving Eq. 8 provide a random complement to the 
Wiener-filtered map such that their sum is a noiseless 
signal consistent with the data and the current power 
spectrum. These artificially created fluctuations aver- 
age out after sufficient iterations leaving us with a re- 
construction of the input signal within the area of the 
primary beam, which we show in Figure 8 as the "Fi- 
nal Mean Reconstructed Signal" . For comparison, the 
"Dirty Map", which is F _1 (H R s + ri), is also shown 
in Figure 8 along with the "Input Signal" which is con- 
structed from the input power spectra shown in blue in 
Figure 4. 

5. CONCLUSION 

In this work the extension of Gibbs sampling to inter- 
ferometric observations of polarized signals has been suc- 
cessfully demonstrated. An example of signal reconstruc- 
tion and inference of CMB power spectra from a moder- 
ately large (n p = 128 2 ) mock data set has been provided. 
The validity of our technique in dealing with realistic 
intcrferometric data, including an incomplete uf-plane 
coverage, finite beam size and baseline-dependent noise, 
has been shown. 

A polarization signal cannot always be uniquely de- 
composed into E and B parts on a cut sky. The non- 
uniqueness of the decomposition causes leakage from the 
E-mode into the much weaker B-mode power spectrum. 
Our input signal maps were generated on a flat sky patch 
with periodic boundary conditions. Because the E-B de- 
composition is unique in this domain, the so-called E-B 
coupling problem did not arise in the recovery of the 
signal realization. Since the Gibbs sampler recovers the 
power spectra of the signal realization, we have a good 
agreement between the "input" and the "mean poste- 
rior" spectra in Figure 4. If the signals had been pro- 
duced as patches cut out from full-sky maps, then the 
"realization" and the "mean posterior" power spectra in 
Figure 4 would have been severely contaminated by E- 
B coupling, which can be easily resolved in the flat-sky 
approximation (Bunn 2002). 

Our Gibbs sampling approach is also applicable to 
more realistic cases of interferometic polarimetry Simula- 
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Figure 4. Mean posterior power spectra for each £-bin are shown in black. Dark and light grey indicate la and 2a uncertainties, 
respectively. The binned power spectra of the signal realization are shown in pink. Blue lines are the input CMB power spectra obtained 
by CAMB for a tensor-to-scalar ratio of T/S = 0.01. 
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Figure 5. Correlation matrices of TE, EE and BB power spectra — only the off-diagonal elements are shown. Correlations and anti- 
correlations between nearby power spectrum bins are the results of having a finite beam width and a finite bin size, respectively. Correlations 
are stronger towards lower signal-to-noise values (from TE to BB) and towards higher ^-values. 
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Figure 6. Marginalized posterior joint distributions of EE and TE power spectra for different ^-bins. Samples of power spectra produced 
by Gibbs sampling have non-Gaussian distributions. 
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Figure 7. Wiener filtered maps, a) Temperature, b) Stokes Q and c) Stokes U components of the solution of Eq. 7; (F 'Ri), transformed 
into Stokes variables T, Q and [/ and averaged over all iterations. The Wiener filtered maps provide the information content of the data. 
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Figure 8. Signal maps, a) The signal realization, which is constructed from the input power spectra shown in blue in Figure 4, is 
used as the input map for the interferometer simulation, b) The final mean posterior map is the sum of solutions of Eq. 7 and Eq. 8; 
{F~ 1 R(x + y)), transformed into Stokes variables T, Q and U and averaged over all iterations. It provides the reconstruction of the 
noiseless input signal by the Gibbs sampler within the area of the primary beam, c) The dirty map is simply the inverse Fourier transform 
of the data. The three rows show, from top to bottom, temperature, Stokes Q and Stokes U parameters. 



tions, such as close-packed arrays with systematic errors. 
Such simulations with varying systematics can provide 
an idea about the limitations of interferometers for up- 
coming missions, such as the QUBIC experiment, which 
aims to detect i?-mode polarization anisotropies of the 
CMB signal. 
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